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ABSTRACT 

We consider methods for finding high-precision approximations to simple zeros of 
smooth functions. As an application, we give fast methods for evaluating the el- 
ementary functions log(x), exp(x), sin(2;) etc. to high precision. For example, if x 
is a positive floating-point number with an n-bit fraction, then (under rather weak 
assumptions) an n-bit approximation to log(x) or exp(x) may be computed in time 
asymptotically equal to 13M(n) log2 n as n — )• oo, where M{n) is the time required 
to multiply floating-point numbers with n-bit fractions. Similar results are given for 
the other elementary functions, and some analogies with operations on formal power 
series are mentioned. 

1 Introduction 

When comparing methods for solving nonlinear equations or evaluating functions, it is custom- 
ary to assume that the basic arithmetic operations (addition, multiplication, etc.) are performed 
with some fixed precision. However, an irrational number can only be approximated to arbitrary 
accuracy if the precision is allowed to increase indefinitely. Thus, we shall consider iterative pro- 
cesses using variable precision. Usually the precision will increase as the computation proceeds, 
and the final result will be obtained to high precision. Of course, we could use the same (high) 
precision throughout, but then the computation would take longer than with variable precision, 
and the final result would be no more accurate. 



Assumptions 

For simplicity we assume that a standard multiple-precision floating-point number representation 
is used, with a binary fraction of n bits, where n is large. The exponent length is fixed, or may 
grow as o{n) if necessary. To avoid table-lookup methods, we assume a machine with a finite 
random-access memory and a fixed number of sequential tape units. Formally, the results hold 
for multitape Turing machines. 

^First appeared in Analytic Computational Complexity (edited by J F Traub), Academic Press, New York, 
1975, 151-176. Retyped with minor corrections and postscript by Frances Page at Oxford University Computing 
Laboratory, 1999 (urls updated 2005). 

Copyright © 1975-2010, R. P. Brent. rpb028 typeset using mgX. 
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Precision n Operations 

An operation is performed with precision n if the operands and result are floating-point numbers 
as above (i.e., precision n numbers), and the relative error in the result is 0(2~"). 

Precision n Multiplication 

Let M{n) be the time required to perform precision n multiplication. (Time may be regarded 
as the number of single-precision operations, or the number of bit operations, if desired.) The 
classical method gives M(n) = O(n^), but methods which are faster for large n are known. 
Asymptotically the fastest method known is that of Schonhage and Strassen [71] , which gives 

(1.1) M{n) = 0(n log(n) log log(n)) as n — )■ oo . 

Our results do not depend on the algorithm used for multiplication, provided M{n) satisfies the 
following two conditions. 

(1.2) n = o(M(n)) , i.e., lim n/M{n) = ; 

n— >-oo 

and, for any a > 0, 

(1.3) M(an) ~ aM(n) , i.e., lim ^/""j = 1 . 
^ ^ ^ ^ V y > > aM{n) 

Condition (1.2) enables us to neglect additions, since the time for an addition is 0(n), which is 
asymptotically negligible compared to the time for a multiplication. Condition (1.3) certainly 
holds if 

M(n) ~ cn[log(n)]''[loglog(n)]''' , 

though it does not hold for some implementations of the Schonhage-Strassen method. We need 
(1.3) to estimate the constants in the asymptotic "O" results: if the constants are not required 
then much weaker assumptions suffice, as in Brent [75a,b]. 

The following lemma follows easily from (1.3). 
Lemma 1.1 

If < a < 1, M(n) = for n < 1, and ci < < C2, then 

oo 

ciM(n) < ^M{a^n) < C2M{n) 

k=0 

for all sufficiently large n. 

2 Basic Multiple-precision Operations 

In this section we summarize some results on the time required to perform the multiple-precision 
operations of divison, extraction of square roots, etc. Additional results are given in Brent [75a]. 
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Reciprocals 

Suppose a 7^ is given and we want to evaluate a precision n approximation to 1/a. Applying 
Newton's method to the equation 

f{x) = a - 1/x = 

gives the well-known iteration 
where 

Since the order of convergence is two, only k ~ log2 n iterations are required if xq is a reasonable 
approximation to 1/a, e.g., a single-precision approximation. 

If £fc_i = 0(2~"/^), then Sk = 0(2~"), so at the last iteration it is sufficient to perform the 
multiplication of Xk-i by £k-i using precision n/2, even though axk-i must be evaluated with 
precision n. Thus, the time required for the last iteration is M{n) + M(n/2) + 0{n). The time 
for the next to last iteration is M(n/2) +M(n/4) + 0(n/2), since this iteration need only give an 
approximation accurate to 0(2~"/^), and so on. Thus, using Lemma 1.1, the total time required 
is 

I{n) - (1 + i)(l + 5 + i + • • ■)M{n) - 3M{n) 

as n — >■ GO. 



Division 



Since b/a = b{l/a), precision n division may be done in time 

D{n) ~ 4M(n) 

as n — >■ DO. 



Inverse Square Roots 

Asymptotically the fastest known method for evaluating a~2 to precision n is to use the third- 
order iteration 

where 

£i = axj — 1 . 

At the last iteration it is sufficient to evaluate ax'j to precision n, ej to precision n/3, and 
Xi{£i — fef) to precision 2n/3. Thus, using Lemma 1.1 as above, the total time required is 

Q(n) ~ (2 + i + |)(1 + 1 + 1 + . . .)M(n) ~ |M(n) 

as n — >■ GO. 
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SquEire Roots 

Since 



a.a 2 if a > 



if a = 0, 

we can evaluate \/a to precision n in time 

R{n) ~ ^M(n) 

as n — > oo. Note the direct use of Newton's method in the form 

(2.1) Xi+i = ^{xi + a/xi) 
or 

(2.2) Xi+i = Xi + 



a — X; 

is asymptotically slower, requiring time ~ 8M{n) or ~ 6M(n) respectively. 

3 Variable- precision Zero-finding Methods 

Suppose C 7^ is a simple zero of the nonlinear equation 

fix) = . 

Here, f{x) is a sufficiently smooth function which can be evaluated near (, with absolute error 
0(2~"), in time w(n). We consider some methods for evaluating to precision n. Since we are 
interested in results for very large n, the time required to obtain a good starting approximation 
is neglected. 

Assumptions 

To obtain sharp results we need the following two assumptions, which are similar to (1.2) and 
(1.3): 

(3.1) M(n) = o{w{n)) , i.e., lim. M{n)/w{n) = ; 

n— >oo 

and, for some a> 1 and all P > 0, 

(3.2) w{/3n) ~ /3"w{n) 
as n — >■ DO. 

Prom (3.1), the time required for a multiplication is negligible compared to the time for a function 
evaluation, if n is sufficiently large. (3.2) implies (3.1) if a > 1, and (3.2) certainly holds if, for 
example, 

w{n) ~ cn"[log(n)]'''[loglog(ra)]'' . 
The next lemma follows from our assumptions in much the same way as Lemma 1.1. 
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Lemma 3.1 

If < /3 < 1, w{n) = for n < 1, and 

Ci < 1/(1 - /3") < C2 , 

then 

oo 

ciw{n) < '^^w{P''n) < C2w{n) 
k=o 

for all sufficiently large n. 

A Discrete Newton Method 

To illustrate the ideas of variable-precision zero-finding methods, we describe a simple dis- 
crete Newton method. More efficient methods are described in the next three sections, and in 
Brent [75a]. 

Consider the iteration 

Xi+i = Xi- f{xi)/gi , 
where gi is a one-sided difference approximation to f'{xi), i.e., 

f{xi + hi) - f{xi) 



9i 



hi 



If £ = \xi — (\ is sufficiently small, f{xi) is evaluated with absolute error 0{e^), and hi is small 
enough that 

(3.3) = f{xi) + 0{ei) , 

then the iteration converges to C, with order at least two. To ensure (3.3), take hi of order £j, 
e.g. hi = f{xi). 

To obtain ^ to precision n, we need two evaluations of / with absolute error 0(2~"'), preceded 
by two evaluations with error, 0(2""/^), etc. Thus the time required is 

(3.4) t{n) ~ 2(1 + 2" + 2-2" + • • ■)w{n) . 
Asymptotic Constants 

We say that a zero-finding method has asymptotic constant C{a) if, to find a simple zero C 7^ 
to precision n, the method requires time t{n) C{a)w{n) as n — t- oo. (The asymptotic constant 
as defined here should not be confused with the "asymptotic error constant" as usually defined 
for single-precision zero-finding methods.) 

For example, from (3.4), the discrete Newton method described above has asymptotic constant 

Civ(a) = 2/(1 - 2-) < 4 . 

Note that the time required to evaluate to precision n is only a small multiple of the time 
required to evaluate f{x) with absolute error 0(2^"). (If we used fixed precision, the time to 
evaluate ( would be of order log(n) times the time to evaluate f{x).) 
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4 A Variable- precision Secant Method 



The secant method is known to be more efficient than the discrete Newton method when fixed- 
precision arithmetic is used. The same is true with variable-precision arithmetic, although the 

ratio of efficiencies is no longer constant, but depends on the exponent a in (3.2). Several secant- 
like methods are described in Brent [75a]; here we consider the simplest such method, which is 
also the most efficient if a < 4.5243 . . . 

The secant iteration is 



fi fi 



i-1 



where fi = f{xi), and we assume that the function evaluations are performed with suficient 
accuracy to ensure that the order of convergence is at least p = {l + ^/5)/2 = 1.6180 . . . , the 
larger root of 

(4.1) />' = P + 1. 

Let s = l^i-i — CI- Since the smaller root of (4.1) lies inside the unit circle, we have 
and 

To give order p, fi must be evaluated with absolute error 0(e^ ). Since fi = 0{\xi — (^\) = 
0(£''), it is also necessary to evaluate {fi — fi-i)/{xi — .Tj_i) with relative error 0(e^ -P), but 
\xi — Xi-i\ ~ e, so it is necessary to evaluate with absolute error 0{£'^ ~'^~^^). [ Since fi must 
be evaluated with absolute error 0{eP ), fi-i must be evaluated with absolute error 0{ef), but 
— /9 + 1 = 2 > p, so this condition is superfluous.] 

The conditions mentioned are sufficient to ensure that the order of convergence is at least p. 

2 

Thus, if we replace by 2"", we see that C may be evaluated to precision n if / is evaluated with 
absolute errors 0(2""), 0(2-2"^"'), 0(2-2"^"'), 0(2-2"''"*), j^. foUows that the asymptotic 
constant of the secant method is 

Cs{a) = 1 + (2p-2)«/(l - P-) < Csil) = 3 . 

The following lemma states that the secant method is asymptotically more efficient than the 
discrete Newton method when variable precision is used. 

Lemma 4.1 

Cs{(x) < CAr(a) for all a > 1. In fact Cs{a)/CN{oi) decreases monotonically from | (when 
a = 1) to ^ (as a — )• oo). 

5 Other Variable-precision Interpolatory Methods 

With fixed precision, inverse quadratic interpolation is more efficient than linear interpolation, 
and inverse cubic interpolation is even more efficient, if the combinatory cost (i.e., "overhead") 
is negligible. With variable precision the situation is different. Inverse quadratic interpolation 
is slightly more efficient than the secant method, but inverse cubic interpolation is not more 
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efficient than inverse quadratic interpolation if a < 4.6056 . . . Since the combinatory cost of 
inverse cubic interpolation is considerably higher than that of inverse quadratic interpolation, 
the inverse cubic method appears even worse if combinatory costs are significant. 

Inverse Quadratic Interpolation 

The analysis of variable-precision methods using inverse quadratic interpolation is similar to that 
for the secant method, so wc only state the results. The order p = 1.8392 ... is the positive 
root of = + It is convenient to define a = 1/ p = 0.5436 ... To evaluate C, to precision 
n requires evaluations of / to (absolute) precision n, (1 — a + (T^)ra, and a^{l — a — + 2a^)n 
for j = 0, 1, 2, . . . Thus, the asymptotic constant is 

CQ{a) = l + {l-a + a^r + i3a^r/{l-a^) 
< Cq(1) = 1 (7 _ 2(7 - (7^) = 2.8085 .... 

Lemma 5.1 

Cqia) < Cs{a) for all a > 1. In fact, CQ{a)/Cs{a) increases monotonically from 0.9361. . . 
(when a = 1) to 1 (as Q — )■ oo). 

Inverse Cubic Interpolation, etc 

U H = 0.5187 ... is the positive root of /x^ + /x^ + /x^ + /x = 1, then the variable-precision method 
of order l//x = 1.9275 . . . , using inverse cubic interpolation, has asymptotic constant 

Cc-(a) = 1 + (1-m + m2)" + (1-m-M^ + V)" + (V)"/(1-M") 
< Cc(l) = (13 -6^-4/^2 -2^3^/3^2.8438... . 

Note that Cc(l) > Cq{1). Variable-precision methods using inverse interpolation of arbitrary 
degree are described in Brent [75a] . Some of these methods are slightly more efficient than the 
inverse quadratic interpolation method if a is large, but inverse quadratic interpolation is the 
most efficient method known for a < 4.6056 ... In practice a is usually 1, 1^ or 2. 

An Open Question 

Is there a method with asymptotic constant C{a) such that C(l) < Cq{1)7 

6 Variable-precision Methods using Derivatives 

In Sections 3 to 5 we considered methods for solving the nonlinear equation f{x) = 0, using 
only evaluations of /. Sometimes it is easy to evaluate f'{x),f"{x),... once f{x) has been 
evaluated, and the following theorem shows that it is possible to take advantage of this. For an 
application, see Section 10. 

Theorem 6.1 

If the time to evaluate f{x) with an absolute error 0(2""') is w{n), where w{n) satisfies conditions 
(3.1) and (3.2), and (for k = 1,2,...) the time to evaluate F^'^\x) with absolute error 0(2~") 
is Wfe(n), where 

Wk{n) = o{w{n)) 
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as n — >■ oo, then the time to evaluate a simple zero ^ 7^ of f{x) to precision n is 



t{n) ~ w{n) 



as n — >■ 00. 



Proof 



For fixed k > 1, we may use a direct or inverse Taylor series method of order k + 1. The 
combinatory cost is of order A;log(fc + l)M(n) (see Brent and Kung [75]). Prom (3.1), this is 
o{w{n)) as n ^ 00. Thus, 

t{n) < [1- {k + iy]-^w{n) +o{w{n)) 
< (l + i + o(l))«^(n). 

For sufficiently large n, the "o(l)" term is less than 1/k, so 

t{n) < (1 + l)w{n) . 

Given e > 0, choose k >2/e. Then, for all sufficiently large n, 

w{n) < t{n) < (1 + e)w{n) , 

so t{n) ^ w{n) as n — >■ 00. 
Corollary 6.1 

If the conditions of Theorem 6.1 hold, /: [a, b] — )• /, f'{x) / for x G [a, b], and g is the inverse 
function of /, then the time to evaluate g{y) with absolute error 0(2~"), for y e I, is 



as n — >■ 00. 
Note 

Corollary 6.1 is optimal in the sense that, if Wg{n) ~ cw{n) for some constant c < 1, then 
w{n) ~ cwg{n) by the same argument, so w{n) ~ (?w{n), a contradiction. Hence, c = 1 is 
minimal. 

7 The Arithmetic-geometric Mean Iteration 

Before considering the multiple-precision evaluation of elementary functions, we recall some 
properties of the arithmetic-geometric (A-G) mean iteration of Gauss [1876] . Starting from any 
two positive numbers oq and h^, we may iterate as follows: 



Wg{n) ^ w{n) 



tti+i — 



aj + bj 
2 



arithmetic mean 



and 




geometric mean 



for 1 = 0,1 
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Second-order Convergence 

The A-G mean iteration is of computational interest because it converges very fast. If bi <^ ai, 
then 

2y/bi/ai r—- 

J- \ Oil Oii 

so only about | log2(ao/6o)| iterations are required before ai/bi is of order 1. Once and bi are 
close together the convergence is second order, for \ibi/ai = l — £i then 

Ei+i = 1 - bi+i/ai+i = 1 - 2(1 - £i)5/(2 - Si) = eil^ + . 
Limit of the A— G Mean Iteration 

There is no essential loss of generality in assuming that oq = 1 and 6o = cos (f) for some 0. If 
a = lim Oj = lim 6j, then 

(7.1) c 



2K(4,) ' 

where i^(0) is the complete elliptic integral of the first kind, i.e., 

7r/2 



J V 1 - sin"' 



de 



' ^ 
(A simple proof of (7.1) is given in Melzak [73].) 

Also, if Co = sin^, q+i = — Oj+i (i = 0, 1, . . .), then 

where E((p) is the complete elliptic integral of the second kind, i.e., 

7r/2 







Both (7.1) and (7.2) were known by Gauss. 

Legendre's Identity 

For future use, we note the identity 

(7.3) m)E{cl)') + K{4>')E{cl>) - K{4>)K{ct>') = ^tt , 

where 0+0' = ^vr. (Legendre [11] proved by differentiation that the left side of (7.3) is constant, 
and the constant may be determined by letting 0^0.) 
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8 Fast Multiple-precision Evaluation of tt 



The classical methods for evaluating tt to precision n take time O(n^): see, for example, Shanks 
and Wrench [62]. Several methods which are asymptotically faster than O(n^) are known. For 

example, in Brent [75a] a method which requires time 0(M(n) log^(rz)) is described. Prom the 
bound (1.1) on M(n), this is faster than 0{n^~^^) for any £ > 0. 

Asymptotically the fastest known methods require time 0(M(n) log(n)). One such method 
is sketched in Beeler et al [72]. The method given here is faster, and does not require the 
preliminary computation of e. 

The Gauss-Legendre Method 

Taking (p = (p' = 7r/4 in (7.3), and dividing both sides by tt^, we obtain 

(8.1) [2K{^/4)E{7r/A) - K^iTr/A)]/^' = ^ . 

However, from the A G mean iteration with ao = 1 and bo = 2~2, and the relations (7.1) and 

(7.2) , we can evaluate i^(7r/4)/7r and i!J(7r/4)/7r, and thus the left side of (8.1). A division then 
gives TT. (The idea of using (7.3) in this way occurred independently to Salamin [75] and Brent 
[75b].) After a little simplification, we obtain the following algorithm (written in pseudo- Algol): 

A^l; B T ^ 1/4; X ^ 1; 

while A- B> 2 " do 

begin Y ^ A; A^^{A + By, B^ VBY ; 
T - X{A-Y)'^ ; 
X ^2X 

end; 

return A^/T [or, better, (A + 5)V(4r)] . 

The rate of convergence is illustrated in Table 8.1. 

Table 8.1: Convergence of the Gauss-Legendre Method 



Iteration 


Ayr - TT 


TT - (A + 5)V(4T) 





8.6 X 10-^ 


2.3 X 10-^ 


1 


4.6 X 10-2 


1.0 X 10-3 


2 


8.8 X 10-5 


7.4 X 10-9 


3 


3.1 X 10"^° 


1.8 X 10-^9 


4 


3.7 X 10-2^ 


5.5 X 10-4^ 


5 


5.5 X 10-43 


2.4 X 10-^4 


6 


1.2 X 10"^^ 


2.3 X 10"^^^ 


7 


5.8 X lO"^'^'^ 


1.1 X 10-345 


8 


1.3 X 10-34^ 


1.1 X 10-694 


9 


6.9 X 10-^^^ 


6.1 X 10-1393 



Since the A-G mean iteration converges with order 2, we need ~ log2 n iterations to obtain pre- 
cision n. Each iteration involves one (precision n) square root, one multiplication, one squaring. 
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one multiplication by a power of two, and some additions. Thus from the results of Section 2, 
the time required to evaluate tt is ^M(n) log2 n. 



Comments 

1. Unlike Newton's iteration, the A-G mean iteration is not self-correcting. Thus, we cannot 
start with low precision and increase it, as was possible in Section 2. 

2. Since there are ~ log2 n iterations, we may lose 0(loglog(n)) bits of accuracy through 
accumulation of rounding errors, even though the algorithm is numerically stable. Thus, 
it may be necessary to work with precision n+0(loglog(n)). Prom (1.3), the time required 
is still ^M(n) log2 n. 



9 Multiple-precision Evaluation of log(a;) 

There are several algorithms for evaluating log(x) to precision n in time 0(M(n) log(n)). For 
example, a method based on Landen transformations of incomplete elliptic integrals is described 
in Brent [75b]. The method described here is essentially due to Salamin (see Beeler et al [72]), 
though the basic relation (9.1) was known by Gauss. 

If cos(^) = £2 is small, then 

(9.1) K{(f)) = (1 + 0(e)) log(4£-^) 

Thus, taking ao = l,6o = where y = 4e~2j and applying the A-G mean iteration to 
compute a = lima^, gives 

iog(y) = ^(i + 0(y-2)) 

for large y. Thus, so long as y > 2"/^, we can evaluate log(y) to precision n. If log(y) = 0{n) 
then ~ 21og2n iterations are required, so the time is ~ 13M(n)log2n, assuming tt is precom- 
puted. 

For example, to find log(10®) we start the A-G mean iteration with ao = 1 and bo = 4 x 10~^. 
Results of the first seven iterations are given to 10 significant figures in Table 9.1. We find that 
7r/(2a7) = 13.81551056, which is correct. 



Table 9.1: Computation of log(106) 



i 






bi 







1.000000000 




4.000000000 X 10- 


-6 


1 


5.000020000 X 10" 


-1 


2.000000000 X 10" 


-3 


2 


2.510010000 X 10" 


-1 


3.162283985 x 10" 


-2 


3 


1.413119199 X 10" 


-1 


8.909188753 x 10" 


-2 


4 


1.152019037 X 10- 


-1 


1.122040359 x 10- 


-1 


5 


1.137029698 x 10" 


-1 


1.136930893 x 10" 


-1 


6 


1.136980295 x 10" 


-1 


1.136980294 x 10" 


-1 


7 


1.136980295 x 10" 


-1 


1.136980295 x 10" 


-1 



Since log (2) = - Iog(2"'), we can evaluate log(2) to precision n in time ~ 13M(n) log2 n. Suppose 
X G [b,c], where 6 > 1. We may set y = 2^x, evaluate log(?/) as above, and use the identity 
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log(a;) = log(y) - nlog(2) 



to evaluate log(x). Since log(y) 2± nlog(2), approximately log2n significant bits will be lost 
through cancellation, so it is necessary to work with precision n + 0(log(n)). 

If X is very close to 1, we have to be careful in order to obtain log(x) with a small relative error. 
Suppose X = 1 + S. If \S\ < 2~"/'°s(") we may use the power series 

log(l + 6) = 6- 6^2 + 6^3 - ... , 

and it is sufficient to take about log(n) terms. If 6 is larger, we may use the above A-G mean 
method, with working precision n + 0(n/log(n)) to compensate for any cancellation. 

Finally, if < x < 1, we may use log(x) = — log(l/x), where log(l/x) is computed as above. To 
summarize, we have proved: 

Theorem 9.1 

If a; > is a precision n number, then log(x) may be evaluated to precision n in time 
~ 13M(n) log2 n as n ^ oo [ assuming tt and log(2) precomputed to precision n + 0{n/ log(n))]. 

Note 

The time required to compute log(x) by the obvious power series method is 0{nM(n)). Since 
13 log2 n < n for n > 83, the method described here may be useful for moderate n, even if the 
classical 0{n^) multiplication algorithm is used. 

10 Multiple-precision Evaluation of exp(a:) 

Corresponding to Theorem 9.1 we have: 
Theorem 10.1 

If [a, h\ is a fixed interval, and x G [a, h] is a precision n number such that exp(a;) does not 
underflow or overflow, the exp(x) can be evaluated to precision n in time ~ 13M(n)log2n as 
n ^ GO (assuming tt and log(2) are precomputed). 

Proof 

To evaluate exp(x) we need to solve the equation f{y) = 0, where f{y) = \og{y) — x, and x is 
regarded as constant. Since 

= (-if-^ {k-l)\y-^ 

can be evaluated in time 0(M(n)) = o(M(n) log(n)) for any fixed A; > 1, the result follows from 
Theorems 6.1 and 9.1. We remark that the {k + l)-th order method in the proof of Theorem 6.1 
may simply be taken as 

k 

Vi+i = yi^{x- \og{yi)y/j ! 

j=0 
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11 Multiple-precision Operations on Complex Numbers 



Before considering the multiple-precision evaluation of trigonometric functions, we need to state 
some results on multiple-precision operations with complex numbers. We assume that a precision 
n complex number z = x + iy is represented as a pair {x, y) of precision n real numbers. As 
before, a precision n operation is one which gives a result with a relative error 0(2~"). (Now, 
of course, the relative error may be complex, but its absolute value must be 0(2~").) Note that 
the smaller component of a complex result may occasionally have a large relative error, or even 
the wrong sign! 

Complex Multiplication 

Since z = (t + iu){v + iw) = {tv — uw) + i{tw + uv), a complex multiplication may be done 
with four real multiplications and two additions. However, we may use an idea of Karatsuba 
and Ofman [62] to reduce the work required to three real multiplications and some additions: 
evalute tv, uw, and {t + u){y + w), then use 

tw + uv = {t + u){v + w) — {tv + uw) . 

Since |i + ?x| < y/2 \t + iu\ and |v + it;| < \/2 \v + iw\, we have 

\{t + u){v + w)\ < 2\z\ . 

Thus, all rounding errors are of order 2~"'\z\ or less, and the computed product has a relative 
error 0(2"'"). The time for the six additions is asymptotically negligible compared to that 
for the three multiplications, so precision n complex multiplication may be performed in time 
- 3M(n). 

Complex Squares 

Since {v + iw)^ = (v — w){v + w) + 2ivw, a complex square may be evaluated with two real 
multiplications and additions, in time ~ 2M{n). 

Complex Division 

Using complex multiplication as above, and the same division algorithm as in the real case, we 
can perform complex division in time ^ 12M(n). However, it is faster to use the identity 

lii^ = (^;2 + w'^y'^Ut + iu)(v - iw)] , 
v + iw 

reducing the problem to one complex multiplication, four real multiplications, one real reciprocal, 
and some additions. This gives time ~ lOM(n). For complex reciprocals we have t = l,u = 0, 
and time ~ 7M(n). 

Complex Squaire Roots 

Using (2.2) requires, at the last iteration, one precision n complex squaring and one precision 
n/2 complex division. Thus, the time required is ~ 2(2 + 10/2)M(n) = 14M(n). 
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Complex A— G Mean Iteration 

Prom the above results, a complex square root and multiplication may be performed in time 
~ 17M(n). Each iteration transforms two points in the complex plane into two new points, and 
has an interesting geometric interpretation. 

12 Multiple-precision Evaluation of Trigonometric Functions 

Since 

(12.1) log(f + iu;) = log |f + iu)| + i.artan(u)/u) 
and 

(12.2) exp{ie) = cos{9) + i. sm{9) , 

we can evaluate artan(a;), cos(a;) and sin(x) if we can evaluate log(z) and exp(2;) for complex 
arguments z. This may be done just as described above for real z, provide we choose the correct 
value of ^/ojbj. Some care is necessary to avoid excessive cancellation; for example, we should 
use the power scries for sin(,T) if \x\ is very small, as described above for log(l + 5). Since 
^ 21og2ra A-G mean iterations arc required to evaluate log(z), and each iteration requires 
time ~ 17M(n), we can evaluate log(2;) in time ~ 34M(n)log2n. Prom the complex version of 
Theorem 6.1, exp(2:) may also be evaluated in time ~ 34M(n) log2 n. 

As an example, consider the evaluation of log(z) for z = 10^(2 + i). The A-G mean iteration is 
started with oq = 1 and bo = 4/z = 1.6 x 10~^ — (8.0 x 10~^)i. The results of six iterations are 
given, to 8 significant figures, in Table 12.1. 

Table 12.1: Evaluation of log 10^(2 + i) 



J 





















(1.0000000, 




0.0000000) 




(1.6000000 X 10" 


-6 


-8.0000000 X 10" 


-7\ 


1 


(5.0000080 X 10- 


-1 


-4.0000000 X 10- 




(1.3017017 X 10- 


-3 


-3.0729008 X 10- 


-4\ 


2 


(2.5065125 X lO' 


-1 


-1.5384504 X 10' 




(2.5686505 x 10- 


-2 


-2.9907884 x 10" 


-3\ 


3 


(1.3816888 X 10" 


-1 


-1.5723167 X 10" 




(8.0373334 x 10" 


-2 


-4.6881008 X 10" 


-3\ 


4 


(1.0927111 X 10" 


-1 


-3.1302088 X 10" 




(1.0540970 X 10" 


-1 


-3.6719673 x 10" 


-3\ 


5 


(1.0734040 X 10" 


-1 


-3.4010880 X 10- 




(1.0732355 X 10- 


-1 


-3.4064951 X 10- 


-3\ 


6 


(1.0733198 X 10- 


-1 


-3.4037916 X 10' 




(1.0733198 X 10- 


-1 


-3.4037918 X 10" 


-3\ 



We find that 

— = 14.620230 + 0.46364761Z 

~ log 1^1 + i.artan(i) 

as expected. 

Another method for evaluating trigonometric functions in time 0(M(n) log (n)), without using 
the identities (12.1) and (12.2) is described in Brent [75b]. 
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13 Operations on Formal Power Series 

There is an obvious similarity between a multiple-precision number with base (5: 

n 
i=l 

and a formal power series: 

oo 

ajX* (oj real, x an indeterminate) . 

Thus, it is not surprising that algorithms similar to those described in Section 2 may be used to 
perform operations on power series. 

In this section only, M(n) denotes the number of scalar operations required to evaluate the first 
n coefficients cq, . . . , c„_i in the formal product 

(oo \ / oo \ oo 

i=0 / \i=0 / i=0 

Clearly, Cj depends only on oq, • • • , and bo, . . . ,bj, in fact 

j 

Cj = dibj-i . 

i=0 

The classical algorithm gives M(n) = 0{n?), but it is possible to use the fast Fourier transform 
(FFT) to obtain 

M{n) = 0(nlog(n)) 

(see Borodin [73]). 

If wc assume that M(n) satisfies conditions (1.2) and (1-3), then the time bounds given in 
Section 2 for division, square roots, etc. of multiple-precision numbers also apply for the corre- 
sponding operations on power series (where we want the first n terms in the result). For example, 

oo 

if P{x) = ^^ajX* and ao ^ 0, then the first n terms in the expansion of 1/P{x) may be found 

i=0 

with ~ 3M(n) operations as n — > oo. However, some operations, e.g. computing exponentials, 
are much easier for power series than for multiple-precision numbers! 

Evaluation of \og{P{x)) 

If ao > we may want to compute the first n terms in the power series Q{x) = log(P(x)). Since 
Q{x) = log(ao) + log(P(a;)/ao), there is no loss of generality in assuming that ao = 1. Suppose 

oo 

Q{x) = From the relation 

1=0 

(13.1) Q'{x) = P'{x)/P{x) , 
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where the prime denotes formal differentiation with respect to a;, we have 

(13.2) J2 = I H ] '''''' 

i=l \i=l / I \i=0 

The first n terms in the power series for the right side of (13.2) may be evaluated with ~ 4M(n) 
operations, and then wc need only compare coefficients to find b, . . . ,bn-i- (Since ag = 1, we 
know that bo = 0.) Thus, the first n terms in log(P(x)) may be found in ~ 4M(n) operations. 
It is interesting to compare this result with Theorem 9.1. 

Evaluation of exp{P{x)) 

If R{x) = exp{P{x)) then R{x) = exp(ao) exp(P(x) — ao), so there is no loss of generality in 

assuming that ag = 0. Now log{R{x)) — P{x) = 0, and wc may regard this as an equation for 
the unknown power series R{x), and solve it by one of the usual iterative methods. For example, 
Newton's method gives the iteration 

(13.3) Ri+iix) = Ri{x) - Ri{x){log{Ri{x)) - P{x)) . 

If we use the starting approximation R{){x) = 1, then the terms in Rk{x) agree exactly with 
those in R{x) up to (but excluding) the term 0{x'^ ). Thus, using (13.3), we can find the first 
n terms of exp(P(x)) in ~ 9M(n) operations, and it is possible to reduce this to ~ ^M(n) 
operations by using a fourth-order method instead of (13.3). Compare Theorem 10.1. 

Evaluation of P"* 

Suppose we want to evaluate (P(x))"* for some large positive integer m. We can assume that 

ao 7^ 0, for otherwise some power of x may be factored out. Also, since P^ = aQ^{P/ao)'^, we 
can assume that oq = 1. By forming P'^,P'^, P^, . . . , and then the appropriate product given by 
the binary expansion of m, we can find the first n terms of in 0(M(n) log2 m) operations. 
Surprisingly, this is not the best possible result, at least for large m. Prom the identity 

(13.4) P™ = exp(mlog(P)) 

and the above results, we can find the first n terms of P'" in 0(M(n)) operations! (If ao 7^ 1, 
we also need 0(log2 m) operations to evaluate a™.) If the methods described above arc used to 
compute the exponential and logarithm in (13.4), then the number of operations is ~ ■^M(n) 
as n — >■ 00. 

Other operations on power series 

The method used to evaluate log(P(x)) can easily be generalised to give a method for f{P{x)), 
where df{t)/dt is a function of t which may be written in terms of square roots, reciprocals 
etc. For example, with /(t) = artan(t) we have df/dt = 1/(1 + t^), so it is easy to evalu- 
ate artan(P(x)). Using Newton's method we can evaluate the inverse function f(-^\P{x)) if 
f{P{x)) can be evaluated. Generalizations and applications are given in Brent and Kung [75]. 

Some operations on formal power series do not correspond to natural operations on multiple- 
precision numbers. One example, already mentioned above, is formal differentiation. Other 

interesting examples are composition and reversion. The classical composition and reversion 
algorithms, as given in Knuth [69], are O(n^), but much faster algorithms exist: see Brent and 
Kung [75]. 
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Shcirper Results 

Some of the constants can be improved. For example, tt can be computed in ~ 6.25M(n) log2 n 
by the Gauss-Legendre method of Section 8, and the constant 13 in Theorem 9.1 can be replaced 
by 10.5. For more information see the postscript to Brent [75a], available electronically from 
http : //wwwmaths . anu . edu . au/~brent/pub/pub032 . html 
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